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Abstract-  Brain  MR  image  segmentation  takes  an  important  role 
in  research  and  clinical  application.  Statistical  method  is 
effective  in  the  segmentation,  which  is  usually  based  on 
maximum  a  posterior  (MAP).  The  key  of  MAP  method  is  to 
estimate  a  prior  probability  of  the  segmentation.  Multilevel 
logistic  (MLL)  model  has  been  used  in  practice  for  the 
estimation.  To  farther  improve  the  performance  of  the 
segmentation,  a  weighted  MLL  (WMLL)  model  is  proposed  in 
this  paper.  The  simulated  results  show  that  the  WMLL  model  is 
effective. 
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I.  Introduction 

Brain  MR  image  segmentation  takes  an  important  role  in 
clinical  applications.  It  is  well  known  that  manual 
segmentation  of  such  images  is  not  only  tedious  but  also 
inconsistent.  So  automatic  or  semi-automatic  methods  are 
desirable.  Hitherto  many  methods  have  been  developed, 
which  can  be  categorized  into  three  classes  [1],  which  are 
based  on  region  [2]-[5],  boundary [6]- [8],  and  point[9]-[17] 
respectively.  The  method  based  on  region  is  computationally 
effective  due  to  the  split  and  merge  algorithm  provided  by 
Horowitz  and  Pavlidis  [2],  but  it  has  difficulty  in  getting  a 
unique  result  [3].  Sonka  [4]  segments  images  into  small 
pieces  and  then  applied  genetic  algorithm  to  rearrange  them, 
so  that  they  form  large  regions  consisting  with  prior 
knowledge,  which  is  obtained  from  manually  segmented 
images.  However,  it  seems  not  easy  to  generalize  such 
method  to  other  cases.  Those  methods  based  on  boundary  fail 
to  sufficiently  utilize  all  messages  in  images.  The 
segmentation  method  based  on  point  can  be  found  in  many 
literatures  recently.  Various  techniques  have  been  adopted, 
including  fuzzy  [9],  neural  networks,  genetic  methods 
[10][11],  statistical  [12]- [17]  and  so  on. 

The  statistical  approach  is  paid  much  more  attentions  in 
present  study,  particularly  those  based  on  maximum  a 
posterior  (MAP)  methods.  In  statistical  methods,  the  prior 
probability  of  the  segmentation  is  not  easy  to  estimate.  So 
maximum  likelihood  (ML)  method  is  applied  in  some 
literature  [17],  which  only  considered  p{y\x}  ,  where  y 
represents  intensities  of  an  image  and  x  the  segmentation.  But 
ML  method  is  liable  to  violate  the  piecewise  congruous  of 
tissues  [17].  As  Markov  random  field  (MRF)  and  Gibbs 
distribution  (GD)  equivalence  [12]  was  introduced,  prior 
distribution  of  the  segmentation  can  be  calculated.  Multilevel 
logistic  (MLL)  model  is  a  typical  way  to  do  so.  To  further 
improve  the  performance  of  the  segmentation,  a  weighted 
MLL  (WMLL)  model  is  proposed  in  this  paper. 

II.  Weighted  multilevel  logistics  model 

There  exist  three  main  obstacles  in  segmentation  of  brain 
MR  images  [17]:  the  thermal  or  electronic  noise,  the  intensity 


non-uniformity  of  same  tissue  classes,  and  the  partial  volume 
effects.  Thermal  noise  is  often  assumed  Gaussian,  white, 
additive,  and  tissue  dependent.  Intensity  non-uniformity  of 
same  tissues  is  due  to  biological  variance  in  different 
structures  of  the  same  tissue  and  irregularities  in  imaging 
equipment.  It  is  assumed  slowly  varied  spatially.  Partial 
volume  effects  are  due  to  the  limited  resolution  of  MR  images. 
Thermal  noise  and  identity  non-uniformity  are  studied  in 
present  research,  but  the  partial  volume  effects  are  left 
unconsidered. 

Let  y={yt  ;  i  ^1}  be  an  image  in  Cartesian  space  N  2, 
where  I EN 2  is  the  region  thaty  occurs.  yt  is  the  intensity  of 
the  image  at  site  i.  y  can  be  regard  as  a  realization  of  analog 
random  variables  at  I.  Suppose  there  are  K  different  tissues 
(classes)  in  the  image,  and  each  of  them  is  labeled  by  a 
number  in  A={1,2,  ...K}.  Let  xt=k ,  k^A  indicates  that  site  i 
belongs  to  class  k ,  then  x={xt ;  i  E 1 }  denotes  a  segmentation 
of  y.  x  can  also  be  regarded  as  a  realization  of  discrete  random 
variables  at  /.  The  goal  of  image  segmentation  is  to  find  an 
optimal  or  sub-optimal  x  under  some  principle. 

In  MAP  method,  P{x\y}  is  considered  and  the  optimal  x=x 
that  makes  it  largest  is  regarded  as  the  segmentation  result. 
Usually  by  Bayes’  theorem  it  can  be  written  in  the  following 
form 

!fP  (1) 

Kp{v\x\p\x\ 

So  the  problem  is  transferred  to  finding  maximum  product 
of p{y\x}  and  P{x}. 

p{y\x}  is  the  joint  density  function  of  y={yi,  y2,  •••,  y\i\} 
under  segmentation  x,  where  |/|  is  the  total  number  of  sites  in 
image  y.  yt  is  assumed  disturbed  by  additive,  white,  Gaussian, 
tissue  dependent,  and  space  variant  noise,  so  if  xf=k,  then 

yi=Mi,k+ni,k  (2) 


where  pik  is  the  mean  intensity  of  tissue  k  at  site  /,  and  nik  is 
Gaussian  noise  at  i  corresponding  to  tissue  k ,  whose  density 


function  obeys  N(0,  o  ik).  Since  the  noise  in  the  image  is  a 
space-variant  white  Gaussian  process  and  tissue  dependent, 
they  are  conditionally  independent  [17].  Let  Rk  be  the  set  of 
sites  that  belong  to  class  k ,  then p{y\x}  can  be  expressed  as 
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where  all  of  parameter  pairs  ={(pih  oik);  k£  A}  form  a 
parameter  set  0={  6  •; i^I}. 

MLL  model  is  an  effective  way  to  estimate  P{xj  by  using 
information  of  the  segmentation  x  itself.  However,  it  still  has 
some  disadvantages  due  to  its  assumption  of  homogeneity  of 
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the  images.  For  example,  different  segmentations  may  have 
the  same  prior  probability  according  to  the  model.  In  fact, 
since  white  matter  (WM)  is  always  surrounded  by  gray  matter 
(GM)  in  brain  MR  images,  any  estimation  of  the  prior 
probability  of  segmentations  that  WM  is  outside  of  GM 
should  be  zero,  or  at  least  very  small,  which  is  not  the  case  in 
MLL  model.  This  implies  that  the  assumption  of 
homogeneous  is  not  acceptable  in  brain  MR  image 
segmentations. 

In  practice,  computation  burden  is  so  large  that  images  are 
segmented  site  by  site  [17].  So  for  site  /,  the  probability  of 
Xj=k  is 

p{xt  =  k\yi,x],j*i)=P{xi  =k\yi,xJ,jsT]i\ 
x  p{yi ,  Xj ,  j  e  n,  I  =  k}p{xi  =k) 

=  p{yi  I  *,•  =  k}p{xj ,  j  e  77,  I  X,.  =  k}p{xi  =  k]  ^ 

=  p{yi  I  xt  =  k}p{xt  =  k,  Xj ,  j  e  ij, } 

=  p{yi\xi=k)Pk{xi} 

provided  that  yi  and  xj’s  are  conditionally  independent,  where 
7/  is  the  neighborhood  of  /,  which  does  not  contain  i  itself, 
and  Pk{Xi}  denotes  P{xt=k,  xJf  j  E  Pj}. 

In  MLL  model,  PkfaJ  is  calculated  as 


where  c  is  a  clique,  Vc(x)  is  the  potential  of  clique  c ,  and  Z  is 
the  normalization  constant.  But  real  images  may  not  be  MRF, 
so  estimating  Pk{xt}  using  (5)  may  produce  errors. 

In  our  study,  Pk{xi}  is  weighted  by  the  probability  of  class 
k  at  site  i  in  order  to  get  better  estimation.  Suppose  that  xk{xi} 
be  the  probability  of  class  k  at  site  /,  then  PkfxJ  can  be  written 
as 


(6) 


where  Z  is  the  normalization  constant. 

Let  t  be  classifications  of  sites  in  7/,  then  Vc(x)  is  exactly 
determined  by  xt  and  t ,  since  in  (5)  or  (6),  i  Ec.  So  Vc(x)  is 
just  Vc(xh  t).  Thus  the  normalization  constant  can  be 
expressed  as 
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where  Q  is  the  space  of  t .  Since  normalization  constant  does 
not  affect  classification,  it  is  denoted  by  Z  for  simplicity  in  the 
context,  despite  that  it  may  have  different  expressions.  So 
P{xt=k  |  yu  Xj,  j  P^i}  can  be  written  as 
p{x;  =k|  yi,Xj,j^i}ocp{yi  I  x;  =  k}Pk{x;} 
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Usually  ft  k{Xi}  is  not  known  except  that  many  accurately 
segmented  images  are  statistically  analyzed.  To  overcome  this 
obstacle,  intensities  of  the  image  are  involved  in  estimating 
7i  k{xi}.  Since  class  k  has  a  mean  intensity  n  k)i  at  site  /, 
k=l,2, ... K ,  it  is  natural  to  assume  that  the  closer  yt  is  to  n  k>i, 
the  more  possible  that  it  belongs  to  class  k ,  and  the  larger 
probability  of  class  k  at  site  i.  An  instinct  way  of  describing 
“closeness”  is  to  adopt  intensity  Euclidean  distance  measure, 
which  is  the  absolute  value  of  difference  between  two 
intensities.  Thus,  the  larger  the  intensity  Euclidean  distance 
between  j/*  and  P  k,b  the  smaller  probability  of  7ik{Xi}.  So 
1  1  (10) 
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Images  are  usually  corrupted  by  thermal  noise,  as 
described  before,  intensity  of  a  single  site  may  be  distorted 
too  much  that  it  cannot  correctly  reflect  intensity  Euclidean 
distance.  So  a  misclassification  may  occur.  To  solve  this 
problem,  two  methods  can  be  applied.  The  first  one  is  to 
decrease  effect  of  noise  in  calculating  7Tk{xJ\  the  second  one 
is  to  add  threshold  in  the  judgment  step  of  the  algorithm. 

To  decrease  effect  of  noise  in  calculating  7ik{xi},  \yr  u  k>l\ 
in  (8)  is  replaced  by  y.  -  juki  ,  which  is  the  mean  value  of 


\yj-  P  k,i\J  ^Qi ,  where  Qt  is  a  neighbor  region  of  site  i.  So 
~  i  i  (ii) 
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This  method  holds  provided  that  the  spatial  probability 
distribution  of  each  class  is  a  smooth  surface. 

A  threshold  qt  can  be  applied  to  decrease  occurrence  of 
the  misclassification  either.  Since  classes  are  piecewise 
contiguous,  if  site  i  belongs  to  some  class  k ,  there  should  be 
some  sites  in  its  neighborhood  7  i  that  also  belongs  to  class  k. 
if  not,  site  i  should  not  be  classified  as  class  k.  qt  takes  the 
role  of  checking  whether  there  is  enough  sites  that  belong  to 
class  k. 

The  above  estimation  does  not  include  noise  model  in 
images,  at  least  it  does  not  care  standard  variance.  In  fact, 
since  noise  is  assumed  to  be  white  and  Gaussian,  this  model 
can  be  applied  in  estimating  probability  of  class  k  at  site  i.  In 
this  case,  the  probability  of  class  k  at  site  i  has  the  following 
form: 
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For  the  same  reason  shown  above,  to  prevent  effect  of 
large  noise  at  some  sites,  \yr  p  k>i\  in  (12)  can  be  replaced  by 


|_y.  _  juk  ,| ,  or  threshold  can  be  introduced,  as  is  done  when  only 
intensity  Euclidean  distance  is  applied. 
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Also  (13)  holds  provided  that  the  spatial  probability 
distribution  of  each  class  is  a  smooth  surface 


III.  Segmentation  algorithm 


Iterated  conditional  modes  (ICM)  proposed  by  Beseg  [13] 
is  applied  to  segment  images  [17].  This  procedure  is  also 
applied  in  our  study.  The  key  point  of  ICM  in  the  study  is  to 
estimate  parameters  in  MLL  model  from  current 
segmentation  result  [14],  then  classifying  the  image  site  by 
site,  then  use  the  new  segmentation  result  to  estimate  MLL 
model  parameters,  and  so  on.  This  procedure  continues  until 
some  criterion  is  satisfied,  say,  maximum  loop  time  is 
reached  or  the  change  between  two  recent  segmentations  is 
ignorable. 

Before  Uk(xi)  is  calculated,  parameter  set  0  should  be 
known.  It  can  be  obtained  by  taking  sample  means  and 
deviations  of  intensities  of  those  sites  belonging  to  class  k  in 
neighbor  region  Qt  of  site  i.  In  case  that  there  are  not  enough 
sites  belonging  to  some  class  k  in  Qh  class  k  is  regarded  as  not 
existed  at  site  i.  So  a  threshold  qt  is  set  so  that  any  class  with 
less  than  qt  sites  belonging  to  it  in  Qt  cannot  be  assigned  to  xt. 

IV.  Results 


Brain  MR  images  are  downloaded  from  brainweb  [18]. 
The  size  of  McGill  data  is  181*217*181  with  voxel  size  of 
1*1*1  mm3,  from  which  the  94th  slice  is  randomly  selected 
for  segmentation.  McGill  provides  two  sets  of  images  with 
different  noise  disturbance  and  intensity  non-uniformity.  One 
is  corrupted  by  3%  noise  with  20%  intensity  non-uniformity, 
and  the  other  is  corrupted  by  9%  noise  with  40%  intensity 
non-uniformity.  The  later  is  adopted  to  check  the  WMLL 
model. 

Since  the  destination  of  segmentation  in  this  study  is  to 
separate  GM  and  WM,  tissues  like  skin,  skull  and  others  are 
roughly  pre-eliminated  by  a  template  before  segmentation. 
According  to  intensity  histogram,  four  kinds  of  tissues  exist  in 
the  remaining  image.  To  speed  up  the  segmentation  process 
and  avoid  background  interference  in  the  process,  sites 
outside  of  the  template  are  kept  from  the  process.  Neighbor 
region  Qt  of  site  i  is  41*41,  and  Qt  is  also  41*41  although 
they  are  not  necessarily  be  the  same.  qt  =  8  and  qt  =  0,1,2 
respectively. 

Fig.l.  shows  the  results  of  the  segmentation  using 
different  n  k(xj)s  described  in  section  II.  (a)  is  the  original 
corrupted  image  and  (b)  is  the  intensity  histogram  of  (a),  (c)- 
(e)  show  the  results  when  (10)  is  adopted  with  qt  =  0,1,2 
respectively;  (f)  shows  result  when  (11)  is  adopted;  (g)-(i) 
show  results  when  (12)  is  adopted  with  qt  =  0,1,2 
respectively;  and  (j)  shows  result  when  (13)  is  adopted.  The 


corresponding  covariance  coefficients  of  WM  and  GM  are 
given  in  Table  1. 

IV.  Discussion  and  conclusion 

In  statistical  segmentation  methods,  MLL  model  is  not 
fully  acceptable  in  brain  MR  image  segmentation.  In  this 


i  j 


Fig.l.  Results  when  different  expression  of  xk  (xj  is  adopted,  (a)  Original 
corrupted  image,  (b)  Intensity  histogram  of  (a),  (c)-(e)  (1 1)  is  adopted  with  qt 
being  0,1,2  respectively,  (f)  (12)  is  adopted,  (g)-(i)  (13)  is  adopted  with  qt 
being  0,1,2  respectively.  G)  (14)  is  adopted. 


TABLE  I 

COVARIANCE  COEFFICIENTS  OF  GM  AND  WM 


^ k(Xf) 

(10) 

(ii) 

(12) 

(13) 

Qi 

0 

1 

2 

0 

0 

1 

2 

0 

GM 

0.9037 

0.9057 

0.9133 

0.9186 

0.8805 

0.8920 

0.9111 

0.9198 

WM 

0.9289 

0.9287 

0.9312 

0.9395 

0.9189 

0.9223 

0.9269 

0.9390 

paper,  WMLL  is  used  to  further  improve  the  performance  of 
the  segmentation.  Two  methods  of  weight  estimation  are 
provided.  (11)  and  (13)  hold  provided  that  the  spatial 
probability  distribution  of  each  tissue  is  smooth.  Since  the 
segmentation  results  are  better  than  others  when  (11)  and  (13) 
are  adopted,  it  seems  that  the  assumption  is  correct  according 
to  our  limited  experiments.  If  only  \yr  u  k>i\  is  adopted  without 
taking  mean  value,  as  in  (10)  and  (12),  setting  qt  >  0  is 
important  to  avoid  small  holes  in  the  results. 
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